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Abstract 
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1 Introduction 



The measurement of the total cross section for electron positron annihilation into hadrons 
allows for a unique test of pertubative QCD. The decay rate T(Z — *> hadrons) provides one 
of the most precise determinations of the strong coupling constant a s . In the high energy 
limit the quark masses can often be neglected. In this approximation QCD corrections to 
R = cr(e + e - — > hadrons) /a(e + e~ — > /i + /i~) have been calculated up to order |TJ, |3|. 
For precision measurements the dominant mass corrections must be included through an 
expansion in m 2 /s. Terms up to order a z s m 2 / s [Q and a^m^/s 2 [H are available at present, 
providing an acceptable approximation from the high energy region down to intermediate 
energy values. For a number of measurements, however, the information on the complete 
mass dependence is desirable. This includes charm and bottom meson production above 
the resonance region, say above 4.5 GeV and 12 GeV, respectively, and of course top 
quark production at a future electron positron collider. 

To order a s this calculation was performed by Kallen and Sabry in the context of QED 
a long time ago |J . With measurements of ever increasing precision, predictions in order 
a 2 are needed for a reliable comparison between theory and experiment. Furthermore, 
when one tries to apply the 0(a) result to QCD, with its running coupling constant, the 
choice of scale becomes important. In fact, the distinction between the two intrinsically 
different scales, the relative momentum versus the center of mass energy, is crucial for a 
stable numerical prediction. This information can be obtained from a full calculation to 
order a 2 only. Such a calculation then allows to predict the cross section in the complete 
energy region where perturbative QCD can be applied — from close to threshold up to 
high energies. 

The two-loop result had been calculated in analytical form M in terms of polyloga- 
rithms. The strategy employed in || was based on a direct evaluation of the imaginary 
part of the vacuum polarization (and thus of the cross section) in a first step followed 
by the calculation of the real part via dispersion relations. In fact, this strategy was also 
applied in the two-loop calculation of the polarization function II for the unequal mass 
case which enters e.g. the p parameter @, || 0. 

The same holds true for a special class of three-loop contributions [ 10 1 to II which are 
present in QED and QCD as well — those induced by massless fermion loop insertions 
into the internal gluon propagator. The imaginary part was calculated in terms of poly- 
logarithms with a method fairly similar to the one employed in order a, with infrared 
cutoffs introduced to calculate separately real and virtual contributions. However, the 
calculation techniques of |TTJ were tailored to the "double bubble" topology and cannot 
be generalized to those which arise for purely photonic or gluonic corrections. 

A complete different approach has been employed for example in the original calcu- 
lations of the 0(a 2 ) corrections to R for massless quarks: The full analytical function 
II(g 2 ) is calculated first, and the imaginary part is obtained as a trivial by-product. The 
pitfalls from infrared divergencies due to the separation into real and virtual radiation 
are thus circumvented. However, the Feynman integrals to be solved are generally not 
available in closed analytical form once the quark masses are different from zero. We will 
demonstrate in this paper that this drawback can be circumvented by numerical methods 
which allow to calculate real and imaginary part simultaneously. Those methods exploit 
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again heavily the analyticity of II(g 2 ). Important ingredients are the leading terms of the 
expansion in the high energy region (—q 2 )/m 2 ^> 1, the Taylor series around q 2 = which 
has been evaluated up to terms of order (q 2 ) 7 and information about the structure of II 
in the threshold region. The method can be tested in the case of the "double bubble" 
diagrams against the known analytical result with highly satisfactory results. Increasing 
the number of terms in the large and small q 2 region, the n(g 2 ) can be reconstructed with 
arbitrary precision. However, the results presented below should be sufficiently precise for 
comparison with experimental results in the forseable future. 

The contributions ~ C F , ~ CaCf and ~ CpTni have to be treated separately since 
they differ significantly in their singularity structure. For each of the three functions 
an interpolation is constructed which incorporates all data and is based on conformal 



mapping and Pade approximation suggested in [12, 13, 14, 15]. Since the result for CpTn, 



is available in closed form the approximation method can be tested and shown to give 
excellent result for this case. Reliable predictions for R to order a 2 and arbitrary m 2 /s 
are thus available. 

In this paper only results without renormalization group improvement and resum- 
mation of the Coulomb singularities from higher orders are presented. Resummation of 
leading higher order terms, phenomenological applications and a more detailed discussion 
of our methods will be presented elsewhere. 

This work decribes the used method in detail and extends the analysis from [TT[ , where 
mainly results were presented. It is organized as follows: The notation and a brief layout 
of our methods are presented in Section |^. The large q 2 behaviour of II is listed in Section 
H|, the ansatz for the threshold singularities in Section |j. Most of the calculational efforts 
were spent on the evaluation of the seven lowest Taylor coefficients of II (q 2 ) presented 
in Section |5]. As a simple application of the result for the first Taylor coefficient the 
contribution of a massive quark to the relation between the QED coupling constant in the 
MS and the on-shell scheme is presented in Section || An approximation which makes use 
of the analyticity structure of II and the results of Section |3] - |5] is constructed in Section 
|7|. The numerical results are presented in Section |8|, together with compact representations 
in terms of simple functions. Section |9] contains a brief summary and conclusions. 



2 Notation and Layout of the Calculation 



In the present approach, originally suggested in |15[ , one calculates real and imaginary 
parts of the vacuum polarization II (q 2 ) and exploits heavily its analyticity properties. The 
physical observable R(s) is related to II(g 2 ) by 

R(s) = 12vrlmn(g 2 + ie). (1) 

It is convenient to define 
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The two-loop coefficient f3\ was calculated in [16 



The number of fermions is denoted by Nf and has to be distinguished from the number 
of light (= massless) fermions ni = Nj — 1. The pole mass of the heavy fermion m is related 
to the MS mass m(/i 2 ) by 
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with L = In /i 2 /m 2 . To describe the singularity structure of II in the region close to 
threshold the perturbative QCD potential Wi 

a v (q 2 ) 



V QC D(q 
a v (q 2 




11 4 

— Ca Trii 

3 3, 



-In 



2 + 3 



-Ca 



(9) 



will become important, where the C^Cp and CpTni contributions have been displayed 
separately. To transform the results from QCD to QED the proper group theoretical 
coefficients Cf = 1, Ca = and T = 1 have to be used. 

In this paper we are only concerned with contributions to II(g 2 ) and R(s) which 
originate from diagrams where the electromagnetic current couples to the massive quark. 
Diagrams where the electromagnetic current couples to a massless quark and the massive 
quark is produced through a virtual gluon have been calculated in |Tj| and will not 
be discussed here. In order a s and a 2 all these amplitudes are proportional to Q 2 , the 
square of the charge of the massive quark. The contribution from the diagrams, depicted 
in Figure p], are of order a s , proportional to and have been calculated analytically 
ISIi: 
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For the a 2 calculation the diagrams can be decomposed according to the colour struc- 
ture. The contributions from diagrams with ni light or one massive internal fermion loop 
will be denoted by CpTriiYii and C F TI1 F with the group theoretical coefficients factored 
out (Figure |]). Purely gluonic corrections (Figure |3]) are proportional to C F or C f Ca- 
The former are the only contributions in an Abelian theory, the latter are characteristic 
for the non-Abelian aspects of QCD. It will be important in the subsequent discussion to 
treat these two classes separately, since they exhibit strikingly different behaviour close 
to threshold. The following decomposition of n(g 2 ) is therefore adopted throughout the 
paper 

n = n(°) + ^Mlcvnw 

7T 

2\\ 2 _ 

2 tt( 2 ) _l n ./^_rr( 2 ) _l tv,.tt( 2 ) _l /^_ttt( 2 ) 




C 2 F U^' + C A C F U^ A + C F T ni Tl\ A > + C F TTlf\ . (11) 

All steps described below will be performed separately for the first three iT 2 ). In fact, 
new information will only be obtained for II^ and 11^. The imaginary part of Imllj 2 ^ 



is known analytically already [HJ. Im!!^ can easily be calculated numerically for the 



following reasons: The amplitude with a massive internal fermion exhibits a two particle 
cut with threshold at 2m which can be calculated analytically ]10|] and a four particle cut 



with threshold at 4m which is given in terms of a two dimensional integral [T(| which can 
be solved easily numerically. This contribution is therefore known already and will not be 
treated in this paper. 

For convenience of the reader we give in the remaining part of this introduction a short 
description of the individual steps of the calculation. 

A crucial input for our procedure is the behaviour of II(g 2 ) for q 2 — > — oo. In Section [5] 
we list the terms of order (m 2 /q 2 ) and (m 2 /q 2 ) 1 for n^ 2 -* in the MS and on-shell scheme. 
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Figure 2: Three-loop diagrams contributing to 11} (inner quark massless) and U F (both 
quarks have the same mass m). 




Figure 3: Purely gluonic contribution to 0(a 2 s ). Diagrams with ghost loops are not 
depicted. 

A very important test for the correctness of the result is provided by the m 4 /s 2 terms 
which are only known for R(s), i.e. the imaginary part of II. 

In Section f| we discuss the behaviour of II(g 2 ) at threshold. For the QED-like cor- 
rections the leading 1/v singularities are resummed via Sommerfeld's formula: R" = 
3vx s /2(l — e~ Xs ) with x s = Cfkoi s /v. An expansion up to 0(a 2 ) delivers the leading 
singularity for Il^g 2 — > Am 2 ) and in combination with the "hard" one-loop vertex cor- 
rection also the subleading logarithmic contribution. With the interpretation of a s in R v 
as ay(4p 2 ), with p being is the momentum of one quark, Eq. @ reproduces on one side 

(2) 

the known threshold term for R t and leads on the other side to an expression for the 

(2) 

leading logarithm of R N a- 

Important information is contained in the Taylor series of IT(g 2 ) around zero. The first 
seven nontrivial terms are computed. The calculation is based on the evaluation of three- 
loop tadpole integrals which is performed with the program package MATAD written in 
FORM |19] . In Section § we describe the method in more detail and present our results. 
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Section |6] is devoted to a discussion of the scheme dependence of the polarization 
function as well as to a related question of the connection between the QED coupling 
constant a in the on-shell and the MS scheme. We explicitly derive the connection up to 
and including the terms of order a 2 a 2 s using the results of the previous section. 

Section [7] contains a description of the procedure which combines the results from the 
different kinematical regions. In a first step a function with a nonsingular behaviour at 
high energies is constructed from the known information about II(g 2 ). Also the logarithms 
at threshold are subtracted. The remaining 1/v singularity for IvJ will be removed by 
multiplication with v. The conformal mapping z = 4a; / (1+uj) 2 of the complex g 2 -plane into 
the interior of the unit circle provides a good starting point for the Pade approximation 
method. The different Pade approximants are denoted by [n/m] and are rational functions 
with a polynomial of degree n in the numerator and m in the denominator. 

The results are discussed and explicit formulae and handy approximations for Ra and 
Rna are given in Section |[ Section |] finally contains our conclusions. 



3 Large q 2 Behaviour 



The expansion of II(g 2 ) for large q 2 provides an important constraint on II in the entire 
complex plane. For the imaginary part it has been argued || that the inclusion of the 
constant, plus quadratic plus quartic terms in m provides an adequate approximation 
down to fairly close to threshold. The large q 2 expansion of II(g 2 ) reads as follows (The 
terms proportional to CfT are listed for completeness.): 
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(The logarithmic contributions of the (m 2 /g 2 ) terms were originally obtained in while 
the constants in H"( 2 ) were first published in || §, albeit not in an explicit form 1 . The 
quadratic mass terms were first found in and later confirmed in p2| .) Eqs. (jl4 - 16 ) 
imply the following contributions to R: 
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Note that Rf = Rf> if one considers only m° and m 2 mass corrections. The equality can 
be understood within the renormalization group approach of [|J. 

The quartic contributions to R have been obtained in ||. They will be important for 
comparisons with the present results and read as follows: 
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1 For instance, in Q the relation following Eq. (8) expresses the constant contribution of II^ 2 ) by some 
quantities which are explicitly given in this article. 
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The results have been expressed in terms of the on-shell mass in order to combine 
them easily with the threshold and low q 2 behaviour. As can be seen from Eqs. (0), ( |2"0"| ) 
and (H) the Abelian parts of n^ 2 ) and i? (2) are in this case independent of \x whereas the 
H dependence of non- Abelian and fermionic contributions is connected to the running of 
the coupling constant a s (fi 2 ). For completeness the Eqs. (|T2] - |T7|) are reexpressed in terms 
of the MS mass m: 
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4 Threshold Behaviour 



The constraints on the threshold behaviour of II(g 2 ) originate from our knowledge about 
the non-relativistic Greens function in the presence of a Coulomb potential and its inter- 
play with "hard" vertex corrections. For a theory with non-vanishing (3 function (QED 
with light fermions or QCD) the proper definition of the coupling constant and its running 
must be taken into account. 

For completeness let us recall the 0(a s ) result. To this order R(s) is known analytically 
and reads (expressed in the variable v = yjl — 4m 2 /s): 
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In the last line an expansion for small velocities v is performed. The n 2 /2v term results 
from longitudinal gluon exchange, the potential, whereas the "—4" originates from a 
transversal, hard gluon. 

For the present problem this has to be extended to 0(a 2 ). Let us in a first step discuss 
the Cp terms which are relevant for QCD and QED. For convenience everything will be 
formulated for QED with the trivial translation a — > Cpa s implied for QCD. The leading 
terms in an expansion of x = na/v are given by 
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where B n are the Bernoulli numbers: B\ = 1/6, B 2 = 1/30, B 3 = 1/42, . . .. These terms 
are induced by rescattering through the Coulomb potential. Formally they are obtained 
from the imaginary part of the Greens function. In addition there is a correction factor 
(1 — 4a/ir) from transversal photon (gluon) exchange. In QCD this term is quite familiar 



from the corrections to quarkonium annihilation through a virtual photon p3[ . These 
considerations allow to reconstruct the 1/v and the constant term in order a 2 [15|, |24[ 
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Now it is possible to define the leading term of H A in such a way that its imaginary part 
reproduces Eq. (|38|): 
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In Ref. [15] this information has been included by considering the combination 11^+411^ 
which is free from a logarithmic singularity at q 2 = 4m 2 and we shall adopt the same 
strategy. 

In the discussion of the terms proportional to CfTtii and CaCf additional aspects 
must be taken into consideration if one wants to predict the term of 0(a 2 v °) multiplied 
by powers of logarithms. This is related to the fact that for a prediction of 0(a 2 ) the 
definition of the coupling constant, the renormalization scheme and the choice of scale 
affect the functional dependence of the result. The problem arises already in a discussion 
of QED with additional rii light fermions. Including purely photonic corrections up to 
order a and contributions with ni light fermion loops of order a 2 one obtains in QED 
close to threshold fllO | 
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where the MS definition of the coupling constant has been adopted and the mass of the 
light fermions is set to zero. This result is easily interpreted in the context of the QED 
potential which can be written in the form (q denotes the relative momentum, with q=2p 
in the case at hand): 
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Returning to the MS convention and QCD one now reconstructs the leading terms of 
the real part of Il(g 2 ) close to threshold^] 
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2 The inclusion of the terms proportional to v 2 In v and v In v in addition to those proportional to 
hiv in the ansatz will lead to a significant stabilization and improvement of our prediction. 
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with an unknown constant term cj. This approach can immediately be transferred to 
QCD, replacing in Eq. © a$ ED (q 2 ) by C F af D (Ap 2 ) and using Eq. (g). Thus for the 
CaCf contribution we postulate the following threshold terms: 
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Because R 



(2) 



(2) 

and RjfA develop only a logarithmic singularity for v 



it is not 

possible to determine the leading real part completely, i.e. q and c^a in (HD an d (|46|) 
remain unknown. 



Small qr 2 Behaviour 



The first seven coefficients of the Taylor series around q 2 = are calculated with the help 
of the program M ATAD which has been used previously for the three-loop corrections to 



the p parameter | 2ql , to Ar [2g] and the top loop induced reaction e + e — > |27 



The package is written in FORM [B| . For a given diagram it performs the derivatives with 



respect to the external momentum q to the desired order, raising internal propagators to 
powers up to about 20. It then performs the traces and uses recurrence relations based 



on the integration by parts method [Bq, 2E] to reduce the resulting three-loop tadpole 



diagrams (in dimensional regularisation). We find 
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The coefficients C[ 2) - Cf ] were already presented in ||11|| , C^ 2 \C^' and C^' are new. 
The C 2 term of Cf ) - cf } and cf } are in agreement with []15| and [j3(J], respectively. 



The /i dependence of the terms proportional to CaCf and CpTrii is again due to the 
running of a s (/x 2 ). If the MS definition of the mass is used, the Abelian part also develops 
a n dependence — a consequence of the anomalous mass dimension — and in general 
In 2 /x 2 terms appear. The corresponding expressions can easily be derived using Eq. (|7|) 
and will not be listed. Instead we want to apply the BLM procedure PT| to our results. It 
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4.70 


-1.67 


25.2 


13.5 
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1.65 


5 


0.14 


1.37 


6.57 


4.16 


-1.49 


23.4 


13.0 


0.195 


0.01 


1.12 


6 


0.11 


1.19 


6.43 


3.75 


-1.35 


21.9 


12.5 


0.181 


-0.09 


0.85 


7 


0.09 


1.05 


6.27 


3.41 


-1.24 


20.7 


12.0 


0.169 


-0.15 


0.76 



Table 1: The BLM scale /j,blm is expressed in terms of the original scale m. For the 
numerical values of and Cjp ni — 5 and fi 2 = m 2 , respectively, /i 2 = m 2 is used. The 
double bubble diagrams with two massive fermion loops are also included. C^ LM is the 
coefficient remaining after the (3q term is absorbed. 

suggests to choose the scale of the 0(a s ) term in such a way that the contribution of the 
0(a 2 ) part proportional to (3q is absorbed. This prescription is based on the observation 
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that the remaining coefficients of the a 2 terms are often relatively small. It is possible to 
treat each term of 0{z n ) separately. In Table [I] we list the BLM scale \iblm for each Taylor 
coefficient together with the numerical value of the original (C^) and the 0(a 2 ) term 
which remains after the BLM scale is adjusted (C^ LM ). It is interesting to note that \iblm 
is decreasing with increasing n. This is plausible because the higher Taylor coefficients are 
increasingly dominated by the threshold region and the characteristic scale at threshold 
is the relative three-momentum of the quarks. Note that C^ LM remains nearly constant 
whereas decreases for increasing n (but remember the rapidly decreasing coefficients 

C<°>). 

In comparison to in Table [1] also the corresponding one- and two-loop coefficients 
(C^\C^) and the two- and three-loop coefficients {C^\ C^) in the MS scheme are 
listed. Whereas is roughly of 0(1) varies from approximately 1.9 down to 0.01. 
Also the sign changes. Therefore the BLM procedure is rather unstable and not applicable. 
On the other hand, the MS coefficients C^f 1 are already reasonably small, so there is no 
urgent need for an optimization procedure. 

Using Eq. (j^j it is possible to transform a s (/i 2 ) into av(/x 2 ), which is based on a physical 
definition. We observed that the quantities C„ are already smaller than the corresponding 
on-shell ones (but still much larger than in the MS scheme). After application of the 
BLM prescription, however, the coefficients are larger than the corresponding C^ LM ■ No 
significant reduction of the coefficients is achieved. 



6 The QED Coupling Constant in the MS Scheme 

The vacuum polarization is connected with the physical cross section through the disper- 
sion relation: 

Ylix 1 Jim? s(s — q z ) 16tc z 

As is well known II(g 2 ) by itself is not a physical quantity. This is reflected by the presence 
of an arbitrary subtraction constant Cq in Eq. (ffHf). In QED usually the condition 



11(0) = (49) 
is adopted which obviously leads to Co = 0. Then the photon propagator 

D{q) = 7TT^m (50) 

gets canonically normalized in the vicinity of the photon pole q 2 ~ and the electro- 
magnetic charge e coincides with the one defined in the Thompson limit. On the other 
hand, using the MS scheme to renormalize the vacuum polarization function, the photon 
propagator looks like 

fl <*> = ?I7iW)' (51) 
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which leads to 



a 



a 



1 + e 2 IT(g 2 = 0) 



a 



(52) 



This relation is not well defined if one takes into account massless quark contributions to 
IT (even within perturbation theory), but it correctly describes the part of the relation 
between on-shell and MS electromagnetic coupling constant coming from a massive quark 
coupled to the electromagnetic current. For a sufficiently heavy quark this relation could 
even be used for a reliable theoretical estimation of the quark contribution to the charge 
renormalization. The Taylor series for the vacuum polarization discussed in the previous 
section was evaluated in a first step in the MS scheme, with a constant term given by 
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This leads to the following relation between a and a(fi 2 ): 

3a(/j, 2 ^ 



(53) 



a 



a(/x 2 



a(/i 2 ) 
1 - 



1 - 

3a 



47T 



Cn 



47T 



Co 



(54) 
(55) 



For the QED case (i.e. Cp 
reproduced. 



1,C A = 0,T = 1/2 and N f = 1) the result of |9 



is 



7 Optimized Approximation 

It is evidently a promising strategy to construct, in a first step, a function Tl(q 2 ) whose 
discontinuity across the cut approximates the one of the full vacuum polarization up 
to subleading terms — in other words, which exhibits a similar spectral function. The 
approximation will have different structure for Ua on one side and Una and IT; on the 
other side as a consequence of their different threshold singularities. 
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For the following discussion the on-shell scheme for the mass m is used and the 't Hooft 
For the C\ term we proceed as in [15| and define IT^ as follows: 



mass fi 2 is set to m 2 



n2V) + 4n ( V 

3 



+ 



16vr 2 



1 w 8z J 8z 24 



(56) 



As already mentioned above the combination U^\q 2 ) + 4II^ 1 - ) (g 2 ) is free of logarithmic 
singularities at threshold whereas the 1/v term remains and is treated below. The high 
energy logarithms have to be removed in such a way that the threshold behaviour and 
the polynomial structure for q 2 — > is not destroyed. Therefore we use the function 
G(z) which is essentially the result in the one-loop case. The terms quadratic in G(z) 
remove the l/zln 2 (— z) terms whereas the linear ones develop only single logarithms for 
q 2 — ► — oo. The factor (1—z) ensures that the threshold is not touched. The last two terms 
reestablish the condition that the new polarization function also vanishes for q 2 = 0. 

In a next step we map the complex g 2 -plane into the interior of the unit circle with 
the help of the variable transformation p2) 



1 - J I -q 2 /4m 2 



1 + Jl- q 2 /4m 2 



Am 2 



:i + u;) 2 ' 



(57) 



The special points q 2 = 0,4m 2 , — oo correspond to u = 0, 1, —1, respectively. In Figure |] 
this conformal mapping is visualized. The upper (lower) part of the cut form z = 1 to 
oo is mapped onto the upper (lower) perimeter of the circle. The Taylor series which in 
the variable z only converges for \q 2 \ < 4m 2 is now convergent in the interior of the unit 
circle in the complex to plane. 
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Figure 4: Transformation between the z and uj plane 



To improve the approximation in a way which allows to approach the boundary of the 
domain of analycity it has been suggested to use Pade approximations. Following [ITH] we 
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define 



Pa{u) 



UJ 



ni 2) (g 2 )-n^(-oo 



(2), 



(58) 



The knowledge of n^(g 2 ) in the three different kinematical regions translates into in- 
formation for P(uj) for uj = —1,0 and 1. Due to the fact that P'(— 1) = (this is 
because of the definition of uj) the constant term in the high energy limit, n^ 2) (-oo), 
is subtracted. The factor 1/(1 + uj) 2 then guarantees that in the limit uj — > — 1 the 
power suppressed term of fl^\q 2 ) is projected out whereas the factor (1 — uj) to- 
gether with the l/v-singularity of II(g 2 ) results in a finite value for P(l). For the con- 
struction of the Pade approximants the following information about Pa(oS) is available: 

{Pa(— 1), -Pa(O), P'a{^)i ■■■■> PjPifyi -Pa(I)}- These ten data points are enough to construct 
Pade approximants like [4/5] or [5/4]. 

(2) (2) 

The above procedure has to be modified for II ^ and 11} because of the different 
threshold behaviour. Here the logarithmic singularities of Una/i are subtracted before the 

~ (2) ~ (2) 

Pade approximation is performed. Therefore 11^ and H] looks like 
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Again the function G(z) is used to get rid of the high energy logarithms. The additional 
factor z in front of the (1 — z)(G(z)) 2 part has to be introduced because of the In 2 (-z) 
behaviour of H$ A and U.\ 2 ^ for z - 



NA 

r(2) 



-oo. The terms next to n^'^ r are needed in order 



to ensure that Iv NA ^(0) = 0. The mapping function is defined as 
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(61) 
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and besides Pna/i{~ 1) and Pna/i{0) the first seven derivatives for u = are available. 
Because of the unknown constants q and cna defined in (fi"4l , f46l) P(l) is not available. 
This allows the construction of Pade approximants like [4/4] or [5/3]. 

(2) 

For the non-Abelian part, R NA , there is an other possibility to construct an approxi- 
mation. This second method is based on the observation that for a special choice of the 
gauge parameter £ the complete threshold behaviour and the logarithmic divergency for 
s — > oo is covered by the gluonic double bubble diagrams given by the last two diagrams 
of Fig. H (together with the corresponding ghost graphs) . Recently, exact results for these 
contributions with arbitrary £ were presented p3[ which are repeated in Appendix |B|. It 



is therefore a strong temptation to consider IT^" — n^|^=4 from the very beginning. The 
imaginary part of this difference is already numerically small, vanishes at threshold and 
has no logarithmic singularity for s — > oo. In analogy to Eq. ( |59D and ( |61| ) we define 



P 9 NA^) 
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and proceed in the same way as above. The coefficients of n^(g 2 )|^ = 4 for q 2 
the high energy expansion can be calculated from Il\ 2 \q 2 ) using the formula |£ 
Appendix H) 



«=4 



Hnp) (g 2 )-^nw (g 2 ) 



(62) 
(63) 

-> and 
(see also 

(64) 



in the different kinematical regions. 

We have checked that both methods lead to equivalent results. The second method is, 
however, more appealing and the results presented below are based on the second choice. 



8 Results 



After the Pade approximation is performed the Eqs. (|5q , |61|j63|) are inverted and the thresh- 
old and high energy terms are reintroduced. Thus real and imaginary parts of Il( 2 ) are at 
hand. Let us subsequently discuss the absorptive part . We choose fi 2 = m 2 throughout 
the whole discussion. 

Figure |5] shows the complete result for Ra, Rna and Ri plotted against the velocity 
v (solid line). This presentation expands the threshold region. The dashed lines repre- 
sent the high energy and threshold approximations. The curves are dominated by the 
threshold singularities. In Figure |6| the same results are shown with abscissa x = 2m/y/s. 
This presentation expands the high energy region. In addition to the quadratic approxi- 
mations which are incorporated into R^ also the quartic approximations (Eqs. (p3j - ^5|) ) 



are plotted. They are evidently very well reproduced by our method. The exact result for 
Ri from jnj is indistinguishable from the curves in Figure |5| and |6]. Also different Pade 



approximants [5/4], [4/5], [4/4], .. . deliver practically identical results with a variation 
comparable to the thickness of the lines. Small variations are visible close to threshold 
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after subtraction of the singular and constant term. Figure |7| shows the results of up 
to ten different Pade approximants. The spread indicates the error of the method. The 
agreement in the high energy region is even better. In this Figure only those curves are 
plotted which have no poles inside the unit circle in the u plane. Perfect agreement for 
Ri is observed which is even slightly increased if higher order Pade approximants (solid 
curves) are taken into account. It is hard to detect the exact result which is represented by 
the dotted line. R NA seems to converge to the solid line ([4/4], [5/3] and [3/5]) when more 
moments from small q 2 are included. The dashed lines are from the [3/3], [4/2], [2/4] and 
[3/4], the dotted ones from Pade approximants including only C NA 3 and C NA 4 ([3/1], [2/2] 
and [3/2]). The dash-dotted curve is the [4/3] Pade approximant and has a pole very close 
to uj = 1(1.07 . . .). For the Abelian part a classification of the different results can be seen: 
the dashed lines are [4/2] and [2/4], the solid ones [3/2], [2/3], [5/3], [3/5] and [5/4] Pade 
approximants. 

(2) (2) 

In the following handy approximation formulae for R A {s) and R NA (s) are presented. 
Details can be found in Appendix [A]. The results read: 

„ (2) (1-w 2 ) 4 3tt 4 2619 3 2061 81 / 2 n . 1 - v 




64 8 V J i + v 



+100p 3/2 (l - p) [2.21P (» - 1.57Pi(p) + 0.27P 2 (p)] , (65) 



/?( 2 ) _ d(2) 



351 o297 
+ v V 




«=4 32 32 

+50p 3/2 (l - p) [1.73P (p) - 1.24P x (p) + 0.64P 2 (p)] , (66) 

with p = (1 — v)/(l + v). The first line of Eq. ( |65|) and (|6^ ) consists of the exactly known 
high energy and threshold contribution, whereas the second and third line represents the 
numerically small remainder which is plotted in Figures || and |9|, respectively, together 
with the result from the Pade approximation. 



9 Conclusions and Summary 

Real and imaginary parts of the three-loop polarization of heavy quarks have been calcu- 
lated. These results can be considered as natural extension of results by Kallen and Sabry 
obtained more than four decades ago. At an intermediate step predictions for the moments 
(up to n = 7) are obtained which are basic ingredients for our numerical method, and 
may in addition also be of interest for the evaluation of QCD sum rules involving heavy 
quarks. The results for the imaginary part allow to predict the cross section in order a 2 s . 
The Pade approximation leads to fairly stable predictions with a relative uncertainty for 
the a 2 coefficients of at most 5%. The prediction in its present form can be considered 
as reliable down to fairly small values of v, say to values of x = 4/3ira s /v around 2. For 
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even smaller values of v the leading terms of the expansion in x have to be resummed. 
This aspect will be treated in a future publication. 
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Appendix 

A Approximation Formulae 

(2) (2) 

In this appendix handy approximations for R A (s) and R NA (s) are derived. We write 
R( 2 \s) for the different colour factors in the form: 

R%\ s ) = R^ ana + R^' rem , (A.l) 

J$i00 = R 9 NA A ™ + R 9 £ ),rem - (A.2) 

R^' ana (s) contains the analytically known threshold and high energy behaviour and is 
chosen in such a way that R^' rem (s) vanishes both at threshold and at the high energies. 
In the following the explicit form of _R( 2 )> ana is presented and the numerical approximation 
of i?( 2 )' rem is discussed. 

Let us first consider the C\ part. Taking the imaginary part of Eq. ( |o"6j ) and solving 
formally for R A \s) we arrive at 

rW+ARM = (A.3) 

R A corresponds to the imaginary part of the last term in Eq. fl56|) and contains the 



constant and m? / s terms of the l.h.s. in the limit s — > oo; R A starts at (9(m 4 /s 2 ). Per 
construction R A has a 1/d singularity at threshold which should be subtracted in such 
way that the high energy behaviour is not disturbed. This is achieved with an additional 
factor (1 — t> 2 ) 4 and leads to 

r{ 2) = R M,rern + (l _ v 2 ) i R (2),t h r_ ARil) _m_ 



A A 

j(2),rem 



v 

2619 ,2061 81 / \ . 1 



+ ''—-'"— |(l-^ 2 )lni-^ (A.4) 



3 ' 

64 " 64 ' 8 V" " ) "' 1 + v 



In Figure || R^ A ,Tem (s) is plotted against p = (l — v)/(l+v) for different Pade approxima- 
tions. The task is now to find an analytically simple formula which describes the behaviour 
of this remainder as close as possible. A closer look into the result of the Pade approx- 
imation shows that an expansion of R ( A > ' rem (s) for high energies (p — > 0) has the form: 
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rem (s) = p 3 / 2 (do + dip + d 2 p 2 + ...)• The appearance of the half-integer exponents 
is a relic of the method in particular of the variable u introduced in Eq. (|57D . The p 3 ^ 2 
term imitates the quartic mass correction which also contain In 2 m 2 pieces in an excellent 
way as can be seen from Figure || Thus we performed a fit in the variable p using Leg- 
endre polynomials P n (p) multiplied by a factor p 3//2 (l — p) which ensures the zeros at the 
boundary. In Figure || one can see that for p —>■ a fit with a low order polynomial would 
fail to deliver an acceptable result, so in a first step a suitable (numerical small) term of 
(D(p 3 / 2 ) is added to get a curve which can easily be parametrised including terms up to 
n = 2. The result reads: 



/™2\ 3/2 

R^ rem (s) = -198 f^-j (V-2^ 

+100j9 3/2 (l - p) [2.21P (p) - l-57Pi(p) + 0.27P 2 (p)] , (A.5) 
and is also shown in Figure |8|. 

(2) 

For the non-Abelian part, R^, we exploit that the gluonic double bubble is known 
exactly (see Eqs. (|62|) and Proceeding as above leads to the result: 



R ( n ] a = R 9 £ ] ' rem + RT +v—-v 3 —, (A.6) 

NA NA g £ =4 32 32 ^ ' 



with 



R 9 ^ Tem = -18 (^J (v*-2v 2 y 

+50p 3/2 (l-p) [1.73P (p) - l-24Pi(p) + 0.64P 2 (p)] . (A.7) 



B Exact Result for Double Bubble Contribution 

In this appendix we list for convenience of the reader the exact result of the imaginary 
part of the fermionic and gluonic double bubble diagram. Following |33| R\ , respectively, 
R^P may be written as (x = I, g) 

R%\s) = -i^ln^-P^pW+SP^ 2 ), (B.l) 

with Rq ^ being the moments in the MS scheme with the results: 



R l = -| + ln4, 
12 4^ 16 s V 4 8 s ' 
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is denned in (0), whereas 5^ was originally calculated in [|K| and reads: 

5 m = _(3-^Hi+^) x 

6 

{ Li 3 (p) - 2 Li 3 (l - p) - 3 Li 3 (p 2 ) - 4 Li 3 ( T ^) - 5 Li 3 (l - p 2 ) + H C(3 ) 



+ Li 2 (p) ln( 



1(1 '" K + 2Li 2 (^)ln /l 



2v 2 



- - ln(^^) 36 ln2 lnp - 44 ln 2 p + 49 lnp In (- — ) + In 2 ( 



6 
1 

Inp In v 

2 F 



2C(2) 
1-v 2 



In p — In 



1 -v 

4v 



2 -\ 



36 ln2 + 21 lnp + 16 \nv - 22 ln(l - v 2 ) 

^ { (15 - 6*; 2 - t; 4 ) (Li 2 (p) + Li 2 (p 2 )) + 3 (7 - 22 1; 2 + 7v 4 ) Li 2 (p) 



+ 



+ 



(l-v) (51 -45w-27w 2 + 5w 3 )C(2) 
\-v) (~9 + 33v -9v 2 -I5v 3 + 4v 4 ) 



ln 2 p 



(33 + 22v 2 - 7v 4 ) ln2 - 10 (3 - v 2 ) (1 + v 2 ) hiv 



+ 2 v (3 - v 2 ) In ( 



- (15-22f 2 + 3w 4 ) In 
4 (1 -v 2 



1 - v 
4^2 



In v — 3 In 



lnp 



237 -9Qv + Q2v 2 + 32v 3 -59v\ , 2 . , + v 

H ; lnp — 16 v (3 — v ) In ( 



4 



4 



-2^(39-17^) In(^) 



u (75 - 29 v 2 ) 



(B-3) 
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Figure 5: Complete results plotted against v — y 1 — Am? / s. The high energy approxi- 
mation includes the m 4 / s 2 term. 
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Figure 6: High energy region (x = 2m/\fs). The complete results (full line) are com- 
pared to the high energy approximations including the m 2 / s (dash-dotted) and the m 4 / s 2 
(dashed) terms. 
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Figure 7: Threshold behaviour of the remainder 5R^ for three different Pade approxi- 
mants. (The singular and constant parts around threshold are subtracted.) 
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Figure 8: The remainder R^' rem plotted for different Pade approximants (dashed lines) 



together with the result of the fit (Eq. 



, solid line). 
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